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Abstract. We describe the initial implementation of magnetohydrodynamics (MHD) in our 
astrophysical simulation code GenASiS. Then, we present MHD simulations exploring the 
capacity of the stationary accretion shock instability (SASI) to generate magnetic fields by 
^ adding a weak magnetic field to an initially spherically symmetric fluid configuration that 

—I models a stalled shock in the post-bounce supernova environment. Upon perturbation and 

nonlinear SASI development, shear flows associated with the spiral SASI mode contributes to 
a widespread and turbulent field amplification mechanism. While the SASI may contribute to 
neutron star magnetization, these simulations do not show qualitatively new features in the 
global evolution of the shock as a result of SASI-induced magnetic field amplification. 
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1. Introduction 

Shortly after the discovery of pulsars [M], the potential role of magnetic fields in the core- 
collapse supernova (CCSN) explosion mechanism began to be investigated (e.g., [29 t [3l [36 1 136] ) . 
In principle, a differentially rotating proto-neutron star (PNS) could give rise to magnetically 
powered explosions. An early conclusion, however, was that (unrealistically) rapid rotation and 
C^ strong magnetic fields would be needed at the pre-collapse stage for magnetic fields to play a 

principal role in the explosion dynamics [29\ HH] . 

More recently, interest in strong magnetic fields has returned in connection with a number of 
observables related to CCSNe, including asymmetries in the explosion ejecta [S2], natal neutron 
star kick velocities [28j, and the high-energy electromagnetic activity connected to some neutron 
stars known as magnetars, or Anomalous X-ray Pulsars (AXPs) and Soft Gamma Repeaters 
(SGRs) (e.g., [151 HHI [251 [53] ) . The theoretical discovery of the magneto-rotational instability 
(MRI) and its application to CCSNe [1] has also contributed to renewed interest in the 
possible role of magnetic fields in the explosion of some supernovae (i.e., those from rapidly 
rotating progenitor cores; e.g., [SU [HI [JOl [lOl HJ] ) . Still, magneto-rotationally driven CCSNe 
likely would be peculiar events, since magnetic progenitor cores tend to rotate slowly [23]. 

Leaving aside the explosion mechanism, the relationship between the formation of neutron 
star magnetic fields and CCSNe is still an open and interesting question, particularly in the case 
of AXPs and SGRs [3l] (although magneto-rotational processes are likely involved [38] [6]). 



The lack of sufficient rotational energy in magnetized pre-collapse progenitor cores, as 
predicted by stellar evolution models, has sparked some recent interest in MHD processes in 
non-rotating CCSN environments [121 EOl 113 IZ] • In particular, Endeve et al. |Q21 E] studied 
magnetic field amplification by the stationary accretion shock instability (SASI [Ij). The SASI 
may play an important role in neutrino-powered explosions [H O [371 [35] . Thus, magnetic fields 
may be an important part of a supernova model if the SASI is found to be sensitive to their 
presence. 

This paper complements the investigations initiated in [16] . We begin with a description of the 
numerical methods used (the initial implementation of MHD in our astrophysical simulation code 
GenASiS), including results from some standard tests. Then, with a new set of high-resolution 
simulations, we investigate the growth and impact of magnetic fields during operation of the 
SASI. We also attempt to quantify the levels of neutron star magnetization that may result from 
SASI dynamics. We find that the magnetic energy may grow exponentially in turbulent flows 
driven by the spiral SASI mode. (The magnetic energy growth time, based on the turnover time 
of SASI-driven turbulence, is estimated to be a few milliseconds.) The presence of amplified 
magnetic fields results in less kinetic energy on small spatial scales, but we find no impact of 
magnetic fields on global shock dynamics. However, magnetic field evolution remains sensitive 
to numerical resolution. We argue that MHD processes associated with the SASI may contribute 
significantly to strong, small-scale PNS magnetic fields, and provide a connection between the 
magnetic fields of neutron stars at birth and supernova dynamics. The saturation energies may 
be sufficient to power flaring activity of AXPs, and possibly SGRs. Moreover, their formation 
does not require progenitor rotation. 

2. Equations and numerical methods 

To study magnetic field amplification from SASI-driven flows we solve the equations of ideal 
MHD. The fluid mass, momentum and energy densities obey conservation laws with sources 

f + V^F.S. (1, 

and the evolution of the magnetic field is governed by Faraday's law (the induction equation) 

^ = -VxE, (2) 

where the vector of conserved variables is U = [p, pu, -E] , the vector of fluxes is F = 
[pu, /9uu-|-IP*-BB, (ii^-|-P*)u-B(B-u)]'^, and the vector of sources is S = [0, -pV<l>,/9u-V<&]^. 
The electric field is E = — uxB (assuming a perfectly conducting fiuid). We denote the combined 
vector of evolved variables by W = [U,B]'^. (Thus, F = F(W) and E = E(W).) Variables 
p, w, E = e + pv? /2 + -6^/2, B denote mass density, fluid velocity, fluid energy density, and 
magnetic fleld, respectively. (We adopt units where the vacuum permeability is po = !•) The 
total pressure (thermal plus magnetic) is P* = P + B"^ /2. In this study, the internal energy 
density is related to the pressure via the ideal gas law e = P/(7 — 1), where 7 is the adiabatic 
exponent (ratio of speciflc heats). 

We report on simulations of flows around a compact object (the PNS), and we take the 
gravitational potential to be given by the point-mass formula $ = —GM/r, where G is Newton's 
constant, M the mass of the central object, and r the radial distance from the center. 

2.1. Magnetohydrodynamics scheme in GenASiS 

We adopt an explicit flnite volume method (e.g., [30]) for the simultaneous time- integration of 

Eqs. ([I| and ([2]). Cartesian coordinates are used to represent points in a three-dimensional 



computational domain with sides L^, Ly, and Lz, and volume V. The computational domain is 
subdivided into A^^^^ x Ny x N^ computational cells with sides Ax, Ay, and Az, and volume 
AV. The x-coordinate of the interface separating cells indexed i and i + 1 is denoted 



'1/2 



+ i Ax, with 



0,...,Nx, and the corresponding cell center coordinate is 



Xi = (xj_i/2 + ^j+i/2)/2- Similar definitions apply for the other coordinate dimensions. The 
time step At increments time from t" to t""^"*^ = t" + At. 

Integrating Eq. M over the spacetime volume element AV x At, applying Gauss' theorem, 
and representing time integrals of the fluxes and sources with the rectangle rule results in the 
time-explicit finite volume update formula 
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Angle brackets on U and S (centered on the geometric centers of computational cells) denote 
volume averages, while angle brackets on the fluxes (centered on the faces of computational cells) 
denote area averages. We use the HLL Riemann solver [22] to compute the interface fluxes, given 
left and right states U^ and U^, respectively. The structure between the fastest left and right 
propagating waves (fast magnetosonic waves with propagation speeds denoted —a~ and a^, 
respectively) in the so-called Riemann fan is then represented by a single average state U*. 
The HLL flux (for example in the x-direction F*) is obtained from the Rankine-Hugoniot jump 
conditions cross the waves separating the left and right states: — a~(U^ — U*) = (F^ — F*) and 
a+(U* - U^) = (F*. - F^). The resulting flux is used in Eq. Q; i.e.. 
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where the fluxes f"'^ are evaluated from W"'^*'^'. Similar expressions are obtained for the 
other flux components. The wave speeds are computed from the left and right states of the 
interface [TT] a^ = max[0, ±A^(W^), ibA^(W^)], where A^ = n^; it Cf^x and cj^x is the fast 
magnetosonic speed. The HLL Riemann solver is among the simplest approximate Riemann 
solvers available for MHD. It relies on minimal information about the characteristic structure 
of the underlying hyperbolic system (only two eigenvalues of the flux Jacobian matrix). As 
such, the HLL flux has been frequently used to construct simple, efficient, and robust solvers for 
hyperbolic systems, including classical MHD J33tl54j. and special and general relativistic MHD 
|12| [T9l rni IT3] . The averaging over the Riemann fan results in diffusive evolution of intermediate 
waves (e.g., contact and Alfven). However, the intermediate waves can be restored in the HLL 
framework [50l ESI EH EH EH El] . 

For second-order spatial accuracy the left and right states needed by the Riemann solver 
are reconstructed from cell-centered volume averages with monotonic linear interpolation: 



jjn,L 



(U)^+P,[(U)^](x,+i/2-Xi) and V'^^^^^ = \J^^,-Vx[Vf^^]{x,+i-x,^,/2)- Monotonic 
interpolation of a variable /" is achieved with the the generalized minmod slope (e.g., |26j) 



Vx [fn = Ax-i minmod [^ (/f - f^_,) , 0.5 (/f+i - fli) , ^ (/f+i - /f)] , 



(5) 



where -d G [1, 2]. We use "d = 1.4 in the simulations presented in this paper. 

The magnetic field is evolved with the constrained transport (CT) method of Evans &; Hawley 
|18| . Magnetic field components are collocated on the faces of computational cells. Integrating 
the induction equation (Eq. [2]) over cell faces and At, applying Stoke's theorem, and replacing 



time integrals of electric field components with the rectangle rule results in time-explicit update 
formulae for area averaged magnetic field components 
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Angle brackets on electric field components (centered on the edges of computational cells) denote 
line averages. Note that with the divergence of the magnetic field defined by 



(v-b; 



(B.), 



-yk 



{B.),_ 



ijk 



Ax 



- — \ 



(5.),: 



Ay 



+ 



^^-)^ifc^ 



(^.>, 



ijk- 



Az 



(9) 



the magnetic field update given by Eqs. rt6^-te| results in (V • B)"^ = (V • B)"-^ — independent 
of how the electric fields are computed^ Thus, the magnetic field remains divergence-free, 
provided it is so initially. 

We use HLL-type expressions [33] (see also [27]) to compute edge centered electric fields in 
Eqs. (|6J)-([8J). The 2:-component centered on, say, (xj+i/2i ?/j+i/2) -^fc) is 



{E. 






"'""^ + «x "y -^^ 



n,SE 



"i Oty E. 



n,NE 



{at + ax){ay + a'y 



+ 



Ufj, cx^ 



{at + ax ) 



[Ey 



^n,w] 



(«y + ay ) 



[Bt 



5"'S1 



(10) 



where at and at' are computed by taking the maximum among the four states surrounding 



the edge. The HLL electric field in Eq. (10) considers the four cells sharing the edge where 



the electric field is to be evaluated. It is derived from the principles used in the derivation 
of the HLL fiux in Eq. (ffl), but applied in two dimensions. For the electric field centered 
on {xi^i/2,yj+i/2, Zk), superscripts SW, NW, SE, and NE indicate that a quantity is obtained 
from values in cells with centers {xi,yj,Zk), {xi,yj+i,Zk), {xi+i,yj,Zk), and (xj+i, j/j+i, z^), 
respectively. Magnetic field components are centered on cell faces, and the cells with centers 
(xi, yj, Zk) and (xj, y^+i, z^) share {By)^.^,^^, which is labeled with superscript W. Similarly, the 



cells with centers (xj, yj, z^) and (xj+i, y^-, Zk) share {Bx 



i+l/2jk 



, which is labeled with superscript 



S. Velocity components are cell centered and must be interpolated in two dimensions, while face 
centered magnetic field components are interpolated in one dimension only. The last two terms 



in Eq. ( 10 ) result from averaging over Riemann fans and introduce explicit dissipation of the 
magnetic field. They act to stabilize the time evolution in under-resolved regions, and can 
become dominant if the bulk of the magnetic fields evolve on spatial scales comparable to the 
grid spacing (see Section [3]). 

The update formulae given by Eq. ^ and Eqs. Q-Q are equivalent to the forward Euler 
method and result in first-order temporal accuracy. Second-order accuracy in time is achieved 



Table 1. Li-error and convergence rate for the circularly polarized Alfven wave test. 
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with a total variation diminishing (TVD) Runge-Kutta method (e.g., [44J ) 

(W)'^ = (W)" + £[(W)"] 
(W)"+i = (W)" + (£[(W)"] + £[(W)^)/2, (11) 

which is a convex combination of two forward Euler steps. The spatial operators and sources 
are here represented by £[(W)]. The time step is determined from the Courant-Friedrichs-Lewy 
(CFL) condition: At = C x mm[Atx, ^ty, Atz], where, for example, Atx = Ax/max[a+, a~], 
C < 1 is the Courant number, and the minimum is taken among all computational cells in V. 

We point out that our MHD scheme is similar to the MC-HLL-UCT scheme described in [33j 
and the MHD scheme in the NIRVANA code [54j . (See also the semidiscrete central- upwind 
schemes developed for hyperbolic conservation laws and Hamilton-Jacobi equations in |27j.) 

2.1.1. Numerical tests To demonstrate the correctness of our implementation of the MHD 
solver we present results from two well-known test problems, the circularly polarized Alfven 
wave and the Orszag-Tang vortex, computed in two spatial dimensions. Both tests use periodic 
boundary conditions everywhere. The Courant number is set to C = 0.4, the slope limiting 
parameter is set to i? = 1.4, and adiabatic exponent 7 = 5/3. 

Circularly polarized (CP) Alfven wave This test problem has been used by many authors 
(e.g., [^ [S3l HH]). In particular, Toth [HI] used this test in an extensive comparison study 
of multidimensional MHD schemes. The CP Alfven wave is an analytic nonlinear solution to 
the MHD equations. As such it provides a means to verify the formal order of the scheme 
through convergence testing. The wave propagates with an angle a = tan~^(2) ~ 63.4° relative 
to the X-axis. The problem is set up with constant background density pQ = 1 and pressure 
Pq = O-l- The velocity (magnetic field) components are Ux{Bx) = n||(i?||) cosa — n_|_(i?_|_) sina, 
Uy{By) = ■U||(i?||) sina + u±{B±) cosq, and Uz{Bz) = ^sin(27rx||), with n_|_ = i?_|_ = Asin(27rx|j) 
and xii = xcosQ + ysina. The background magnetic field parallel to the direction of wave 
propagation is B\\ = 1 (i.e., the Alfven speed Ua is unity). The parallel velocity u\\ is set 
to zero for the traveling wave (TW) test and —1 for the standing wave (SW) test. The 
wave amplitude is A = 0.01. Periodic boundary conditions can be used if the domain is 
[x,y\ G [0, 1/cosa] x [0, 1/sina] {Lx/Ly = 2, and we use Nx = 2Ny for square cells). 

Convergence results from the CP Alfven wave test (TW and SW tests), run for one {t = 1) 
and five {t = 5) grid crossings, are listed in Table [11 The Li-error norm is Li{Bz) = Y^ \Bz{t = 
0) — Bz{t = tf)\/Y \Bz{t = 0)1, where the sums extend over all cells. The tests demonstrate 
second-order convergence with increasing resolution. The SW errors are larger (almost a factor 2) 
than the TW errors. (The convergence rate is given by \og{Li{Bz.,Nx)/Li{Bz,Nx/'2))/\og{2).) 



Orszag-Tang Vortex This is another test which has be used extensively to benchmark 

(e.g., 



multidimensional MHD codes (e.g., \3'6\ I45j ) . Following the initialization in 1-^5], the 
computational domain is the unit square [x,y] £ [0, 1] x [0, 1], where the density and pressure 
are initially uniform p = t^ and P = jj^ (resulting in sound speed Cs = 1). The nonzero 
components of the velocity field are Ux = — sin(27ry) and Uy = sin(27ra;), and the magnetic field 
is given by B^ = —Bq sin(27ry) and By = Bq sin(47ra;) with Bq = 1/\/Att. 

Shocks develop quickly in the fiow for i > 0. The flow becomes very complex and develops 
MHD turbulence after multiple shock-shock interactions. Figure [T] displays results from the 
Orszag-Tang test computed with GenASiS. With the proflles displayed in the lower panels, 
we have computed the Li-error norm of the lower-resolution runs with respect to the high- 
resolution run (4096^). The error norm decreases with increasing resolution at the flrst order 
rate (as expected for flows containing discontinuities). 



2.2. Gravitational source terms 

The gravitational force and power must be specified for problems involving gravity. In particular, 
constructing source terms which result in good energy conservation properties is of interest. 
Multiplication of the discrete mass conservation equation in Eq. (|3J) with the potential (^)jjfc 
results in (after some algebra and assuming a static gravitational potential) 
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where the gravitational energy density is {Eg)^-j^ = {p)^jk{^)ijk ^^'^ ^^^ x-component of the 



gravitational energy fiux is (Fx"')r+i/2jfc = {F£)1+i/2jkii^)ijk + {^)i+ijk)/'^- We have also 
defined {FS 9x)'l+i/2jk = iF£)1+i/2jk(i^)i+ijk ~ {^)ijk)/^^- (The x-component of the mass fiux 
in Eq. (p| is denoted {F§)^,-^i2ji^-) The last three terms in Eq. (12) suggest a finite volume 

representation of the gravitational power ("S*^)"-^ = — (/ou • V$)"-^,. This choice results in energy 
conservation to machine precision; i.e., 
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This result relies on the assumption of a static gravitational potential, and is sufficient for the 
simulations presented in this paper. We will present a suitable generalized discretization of the 
gravitational power — valid for time-dependent potentials — in a forthcoming paper. 

A finite volume representation of the x-component of the gravitational force is given by 
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where {P9x)1:^i/2jk = ip)7+i/2jk((^)i+ijk ~ {^)ijk)/^^- (Similar expressions are used for the y- 
and z-components.) For the density centered on x = XJ4.1/2 we use the so-called HLL average 
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Figure 1. Results from running the Orszag-Tang vortex test witli GenASiS. Tlie upper panels 
shows the fluid pressure at times t = 0.5 (left) and t = 1.0 (right) for a model computed 
with 4096^ cells. In the two lower panels we plot the thermal pressure versus a; at i = 0.5 for 
y = 0.4277 (left) and y = 0.3125 (right) (cf. [l5]). Profiles are plotted for runs with varying 
spatial resolutions: A^^ x A^^^ = 128^ ( ), 256^ ( ), 512^ ( ), and 4096^ ( ). 



Note that the expression for the gravitational power does not depend on the particular Riemann 



solver, while the force (with Eq. (15)) does. For a different Riemann solver, Eq. (15) should be 



replaced with the face density consistent with that solver. 



3. Simulations of SASI-induced magnetic field amplification 

In this section we describe high- resolution MHD simulations with GenASiS, exploring the 
capacity of nonlinear SASI-driven flows to amplify magnetic flelds. With an idealized model 
we seek to investigate the character of SASI-induced magnetic fleld ampliflcation. Interesting 
questions to explore include: does the presence of magnetic fields impact the evolution of the 



SASI?, what are the implications for observables associated with CCSNe, in particular PNS 
magnetization? However, with idealized simulations, we can only begin to scratch the surface 
of these complicated topics. 



3.1. Model setup and parallel execution 
Our initial conditions follow closely the adiabatic setup 
described in [4j and [5j (see also ^6]): a spherical 
stationary accretion shock is placed at a radius i?sh = 
200 km. A highly supersonic flow is nearly free-falling 
towards the shock for r > i?sh- Between the shock 
and the PNS the flow settles subsonically — obeying the 
Bernoulh equation u'^/2 + ^P/[p{-f - 1)] - GM/r = 0— 
in nearly hydrostatic equilibrium. Matter is allowed to 
flow through an inner boundary placed at -Rpns = 40 km. 
The ratio of specific heats is set to 7 = 4/3, the mass 
of the central object is held fixed at M = 1.2 M©, and 
the accretion rate ahead of the shock M = 0.36 Mq s^^ 
is also held fixed during the simulations. A radial 
magnetic field is superimposed on the initial condition: 
Br = sign(cos0) xi?o(-RpNs/'")^) where Bq is the magnetic 
field strength at the surface of the PNS. We initiate the 
SASI with random pressure perturbations in the post- 
shock flow. 




Figure 2. Domain composition from 
a 3D MHD-SASI simulation. 
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The MHD solver is parallelized using the Mes- 
sage Passing Interface (MPI), and the computa- 
tional domain is subdivided into blocks containing 
an equal number of zones, which are distributed 
among MPI processes (see Figure ^. During 
the simulations we keep the number of zones per 
block (MPI process) fixed to 32^. To conserve 
computational resources the simulations are initi- 
ated in a relatively small computational domain 
with sides L^-,^-, = 600 km covered by 512 zones 
(Ax ~1.17 km). Once the SASI evolves into the 
nonlinear regime the volume encompassed by the 
shock grows, and the shock eventually interacts 
with the domain boundaries. When this happens, 
we expand the computational domain by adding 
a layer of 32^-zones blocks (i.e., we add 64 zones 
in each coordinate direction) and restart the simu- 
lation from the last checkpoint written before the 
shock interacted with the boundary. We repeat this 
process, and run our simulations until the shock in- 
teracts with the boundary of the largest computational box L^ax = 1500 km, or the simulation 
time reaches t = 1100 ms, whichever occurs first. Since we keep Ax fixed during the simulations, 
the largest computational domain is covered by 1280 zones in each spatial dimension. The weak 
scaling plot in Figure |3] shows the wall-time per time step, which remains reasonably constant as 
the computational domain expands to accommodate the growing shock volume. During a run, 
we write simulation output for analysis and visualization every 2 ms of physical time, resulting 
in tens of Terabytes of data from each model. 
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Figure 3. Weak scaling from a high- 
resolution 3D MHD-SASI simulation. 



3.2. Simulation results 




Figure 4. Volume rendering of 
the Mach number Ma = |u|/cs 
with velocity vectors (|u| > 
10 km s~^) superimposed. An- 
gular momentum redistribution 
has occurred as a result of the 
SASI, and the presence of coun- 
terrotating flows is apparent. 
The shock triple-point [5j (a line 
segment extending across the 
shock surface), positioned to the 
upper right, is visible as the kink 
in the shock surface (white con- 
tour). It connects the pre-shock 
accretion flow and the two coun- 
terrotating post-shock flows, and 
moves on the shock surface in the 
counterclockwise direction. A 
layer of sheared flows extends 
from the triple-point and results 
in post-shock vorticity genera- 
tion. Image courtesy of Ross 
Toedte, ORNL. 



We have carried out simulations with varying initial magnetic field: Bq = 1 x 10^*^ G (model 
BIO or "weak-field model"), Bq = 1 x lO^^ q (model B12), and Bq = I x lO^^ q (model B13 
or "strong- field model"). The initial field can be considered weak in all the models in the sense 
that the magnetic energy density is small compared to the kinetic and internal energy densities. 

Upon perturbation and nonlinear development, all models eventually exhibit typical spiral 
mode fiows. This is consistent with [5], who found the spiral mode to dominate the late-time 
evolution independent of the initial perturbation. Figure |4] illustrates the post-shock fiows 
associated with the nonlinear spiral SASI mode. (The rendering is created from a snapshot of 
the strong-field model at t = 820 ms, but the hydrodynamic developments exhibited by this 
model are typical of all the computed models.) 

The pre-shock accretion fiow impinges on the shock at an oblique angle due to the aspherical 
shock and its off-center position. The sizable tangential velocity component (relative to the 
shock surface), which is preserved across the shock, leads to supersonic post-shock fiows ahead 
of (and directed towards) the triple-point. Supersonic shear fiows are directed down towards the 
PNS and result in vorticity generation. Vorticity is distributed in a large fraction of the post- 
shock volume during the operation of the spiral SASI mode. Strongly forced accretion-driven 
turbulence develops as a result of the SASI, and the post-shock fiow becomes roughly divided 
into a supersonic (driving) component and a subsonic (volume- filling) turbulent component. 
These hydrodynamic developments result in turbulent magnetic field amplification |16l ITT] . 

The magnetic field amplifies in response to the hydrodynamics developments, and the 
magnetic energy becomes concentrated in thin intense magnetic fiux tubes (or ropes). The 
volume rendering in Figure [5] shows streamlines tracing out the magnetic field below the shock 
and reveals a complicated magnetic field geometry. The magnetic field is "frozen-in" to the 




Figure 5. Streamlines tracing magnetic field lines in an MHD-SASI simulation during the 
nonlinear stage. Streamlines are randomly seeded on the shock surface. Image courtesy 
of Dave Pugmire, ORNL. (An animation from which this image is taken can be viewed at 
[http: //events . eels . anl .gov/scidacll/visualization-night/visualization-night -winners/^.) 



fluid in the ideal MHD limit. The large scale flows associated with the SASI result in relatively 
straight field lines, while the turbulent flows are responsible for a more "chaotic" structure. 
Figure p^ contains further details from the numerical simulations. In Panel (A) we plot the 

relative change in total magnetic energy below the shock versus time for all models: BIO ( ), 

B12 ( ), and B13 ( ). (The change is scaled to the initial magnetic energy for easy 

comparison across the models.) After an initial spurt, all the models experience an early period 
of exponential magnetic energy growth with essentially the same growth rate (cf. the temporal 
window from 650 ms to 780 ms). Exponential growth is typical during the kinematic regime of a 
turbulent dynamo, in which the magnetic field's back-reaction on the fluid is negligible (e.g., [7]), 
and the growth time is roughly the turnover time of the turbulent eddies r^diy. The magnetic 
energy in the weak-field model grows exponentially at a nearly constant rate, with growth time 
T ~ 66 ms (dash-dotted reference lines), until the end of the simulation (t = 1100 ms), and 
receives a total boost of about five orders of magnitude {E^ ~ 1.8 x 10^^ B). In the model with 
Bq = Ix 10^^ G (B12), E^ also grows steadily until the end of the run. The magnetic energy in 
this model receives a boost of four orders of magnitude (£'m ~ 2.3 x 10~^ B). It grows initially at 
the same rate as the weak- field model, but the growth rate tapers off at later times {t > 900 ms). 
The strong-field model (B13, Bq = 1 x 10^^ G) exhibits exponential magnetic energy growth 
(r ~ 66 ms) early on. Then, around t ~ 780 ms, the growth rate drops almost discontinuously, 
and Sniag grows by only about 50% for the remainder of the simulation. Model B13 receives a 
total boost in magnetic energy of about a factor of 300 {E^^ ~ 8.9 x 10^^ B). The abrupt drop in 
the magnetic energy growth rate observed in the strong-field model occurs when magnetic fields 



become dynamically important: concentrations of high (3^^i^ ( 



|u| ), which exceeds unity in 



localized regions, are scattered throughout the shock volume. About 50% of the total magnetic 
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Figure 6. Summary of simulation results. Panel (A): magnetic energy versus time; panel (B): 
turbulent kinetic energy versus time; panel (C): growth rates due to various terms in Eq. ( |16| ) 
versus time for model BIO; and panel (D): magnetic rms scale (average flux rope thickness) 
versus time. See text for further details and discussion. (1 B = 10^^ erg.) 



energy in model B13 resides in regions where /3,^„ > 10^^, and these magnetic fields occupy less 
than 10% of the total shock volume. Regions with Wa ^ |u| appear and disappear in a highly 
intermittent manner during this stage. 

The magnetic energy grows at the expense of turbulent kinetic energy below the shock. 
The turbulent kinetic energy E'^ = J"^ e^iWjdk is plotted in Panel (B). It is defined as flows 
varying on scales with wavenumber k > k-j. = 27r/AT, where At ~ 30 km. (The spectral 
kinetic energy density ek(/c) is obtained from Fourier transforms of components of the kinetic 
energy.) The particular choice for k-r is motivated by several factors, including (1) magnetic 



field amplification occurs mostly on spatial scales with k > k^ (from magnetic energy spectra 
we obtain a characteristic magnetic field scale A^ ~ 20 km), and (2) the Taylor microscale 
{{v?)/{\V X up))-*^/^, which measures the average size of turbulent eddies, is comparable to At- 

The turbulent kinetic energy grows rapidly during the ramp-up of the SASI and reaches 
a saturation level. After saturation, the time-averaged turbulent kinetic energy (averaged 
over the time interval from 800 ms to 1100 ms) in the respective models are found to be 
{E^)l-ll = 5.90 X 10-3 B (BIO), 5.48 x lO^^ B (B12), and 5.11 x lO'^ B (B13). (The total 
kinetic energy below the shock E^^ is about an order of magnitude larger during this stage.) 
We note that E"^ in model B12 and model B13 is reduced (relative to the weak-field model) by 
^ 7% (4.2 X 10"^ B) and ~ 13% (7.9 x 10~^ B), respectively. The reduction in turbulent energy 
is comparable to the gain in magnetic energy in the respective models. Although there is a 
measurable reduction in the turbulent kinetic energy in the models attaining stronger S-fields, 
we find no impact of magnetic fields on large scale flows and SASI evolution. (The turbulent rms 
velocity Urms inferred from these models exceeds 4000 km s~^, and the associated eddy turnover 
time is Teddy = >^m/Urms « 5 ms.) 

The saturation level for the turbulent kinetic energy may serve as an upper limit on the 
magnetic energy attainable in these simulations. We note that the magnetic energy in model 
B13 (dotted red line in Panel (B)) seems to saturate at about 10% of the turbulent kinetic energy. 
The magnetic fields in this model become dynamically significant, but are also influenced by 
numerical dissipation during the evolution (see below). The fluid in a developing CCSN is an 
excellent electrical conductor and magnetic energy dissipation is probably not important on the 
explosion time scale. Thus, the magnetic energy may possibly grow beyond the levels seen in 
our simulations, but probably not above the turbulent kinetic energy. 

Assuming a non-ideal electric field — (u x B) + r/J with scalar resistivity rj, the evolution 
equation for the magnetic energy density is easily obtained by dotting Eq. ([2| with B: 

— ^ = B-[(B-V)u-(u-V)B-B(V-u)- V X (r/J)] 

= -V-P-u-(J xB)-B-Vx (r/J), (16) 

where P = [u(B • B) — B(B • u)] and J = V x B. (We use scalar resistivity here for sake of 
simplicity of the discussion, but it is also appropriate due to the similarity between the last two 
terms in Eq. ( |10[ ) and rjJz-) The first and third term on the right-hand-side of Eq. ( |16| ) (first line) 
describe magnetic field amplification due to stretching and compression, respectively. Magnetic 

energy growth rates due to compression ( ) and stretching ( ) are plotted versus time 

for model BIO in Panel (C) of Figure rol Stretching is the dominant magnetic field amplification 
mechanism in the nonlinear regime (i > 650). In a turbulent flow, the separation between 
two (initially adjacent) fluid elements grows exponentially with time. If the fluid elements are 
connected by a weak magnetic field the frozen-in condition of ideal MHD results in stretching, 
and thereby strengthening, of the magnetic field and an increase in the magnetic energy 



The terms on the right-hand-side of Eq. (16) (second line) represent Poynting fiuxes, work 



done against the Lorentz force (VFl) and magnetic energy decay due to resistive dissipation 
{—Qj), respectively. Kinetic energy is converted into magnetic energy if Wi^ > 0. It is apparent 



from Eq. ( 16 ) that the magnetic energy growth rate is due to the sum of individual rates from 
accretion of magnetic energy (Poynting flux) through the surface enclosing the PNS (9VpNs), 
work done against the Lorentz force, and resistive dissipation. The Poynting flux through 9T4>ns 
and resistive dissipation generally result in decay of the magnetic energy in the computational 
domain. The decay must be overcome by the Lorentz work term in order for the magnetic 
energy to increase. (The magnetic Reynolds number is R^ = |Wl|/|Qj|.) Growth rates due to 

work done against the Lorentz force ( ) and losses due to the Poynting flux through 9Vpns 

( — • — ) are plotted in Figure [6] (Panel (C)). Evidently, the growth due to VFl greatly exceeds 



that due to Poynting losses through 91pNs- Moreover, the growth rate is ~ 480 s~^, which 
imphes a magnetic energy growth time of about 2 ms (comparable to r^ddy!)) as opposed to the 
66 ms seen in Panel (A). This discrepancy in growth rates is due to numerical dissipation: the 



field evolves on small spatial scales in the turbulent flows and the dissipative terms in Eq. ( 10 ) 



become nontrivial and affect the overall growth rate. This suggests that the magnetic energy 
may grow on a millisecond time scale in SASI-driven turbulence with high R^. 

The effect of numerical dissipation is illus- 
trated in the evolution of magnetic fiux tubes. 
In Panel (D) of Figure p\ we plot the average 
flux tube thickness A,^, = ((S2)/(|V5|2))V2 
versus time for all the models. During the ramp- 
up of the SASI, Arms decreases rapidly until it 
reaches about 4 km (~ 3 Ax) and stays rela- 
tively constant thereafter. The decrease in flux 
rope thickness is halted by (numerical) resistive 
dissipation when Arms reaches the resistive scale 
(a few X Ax). In an environment with near per- 
fect conductivity (i.e., the CCSN environment) 
we expect the held amplification to be halted by 
dynamical interactions with the fluid before the 
resistive scale is reached |38]. In our large-scale 
simulations, flnite spatial resolution and associ- 
ated numerical dissipation plays a decisive role 
for the magnetic field evolution: the attainable 
field strength and the magnetic energy growth 
rate are underestimated. By investigating indi- 
vidual terms in the induction equation during 
the late-time phase of the MHD-SASI simula- 
tions, we have verified that the dissipative terms 
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in Eq. ( 10 ) result in magnetic energy decay with 



a decay time comparable (but smaller in magnitude) to the growth time due to VFl. 

The geometric concept of magnetic flux ropes is further illustrated in Figure [7j where we plot 
the time-averaged rms magnetic held (-Brms) versus time-averaged magnetic rms scale (Arms) from 
models where the spatial resolution has been varied (Ax ~ 2.34, 1.56, and 1.17 km). Higher 
spatial resolution accommodates the development of thinner flux ropes. The decreased cross- 
sectional area of a flux rope implies an increased field strength. In Figure [7j the increase in field 
strength is in direct proportion to the decrease in the cross-sectional area afforded by higher 
spatial resolution. Also, the slope observed in the 3D simulations is steeper than that observed in 
2D (axisymmetric) simulations [16] . (We also observe that the exponential growth rate increases 
with increasing spatial resolution.) 

The increase in magnetic energy in the volume occupied by the PNS due to Poynting flux 
through aVpNs is 1.1 x 10^^ erg (BIO), 3.2 x 10^^ erg (B12), and 4.5 x 10^^ erg (B13). The 
results form models B12 and B13 imply PNS rms magnetic fields exceeding 10^^ G. We expect 
the degree of PNS magnetization to be insensitive to the initial field in realistic settings with 
high i?m where dissipative processes can be ignored. 



4. Summary and discussion 

We have presented the initial implementation of a second-order finite volume MHD scheme in 
GenASiS. Using the MHD capabilities, we have performed 3D simulations of the evolution 
and amplification of magnetic fields in SASI-driven flows. The simulations are initiated from a 



configuration which resembles the early stalled shock phase in a CCSN, albeit with simplified 
physics that excludes critical components of a supernova model (e.g., neutrino transport, self- 
gravity, and the PNS itself). On the other hand, our simulations are computed with a spatial 
resolution that is currently inaccessible to state-of-the-art supernova models in three spatial 
dimensions, and they may therefore provide valuable insight into MHD developments in CCSNe. 

Flows associated with the spiral SASI mode result in vigorous turbulence below the shock 
("Urms ^ 4000 km s~^). The turbulence amplifies magnetic fields by stretching, and the magnetic 
energy grows exponentially with time as long as the kinematic regime obtains. The resulting 
magnetic fields display a highly intermittent flux rope structure. In the strong-field model, 
the magnetic energy saturates when the associated energy density becomes comparable to the 
kinetic energy density (i.e., Va ~ |u|) in localized regions of the flow. The subsequent magnetic 
field evolution remains highly dynamic: strong fields emerge, are advected with the flow, are 
temporarily weakened, and then reemerge in a seemingly stochastic manner. 

The presence of amplifled magnetic fields does not result in noticeable effects on the global 
shock dynamics in our simulations, and this may be understood from considerations of the 
energetics: the turbulent kinetic energy — which powers SASI-driven field amplification — only 
amounts to about 10% of the total kinetic energy below the shock (i^^ ~ 5 x 10~^ B), and is 
not significant to the explosion (~ 1 B). These observations alone suggest a rather passive role 
of magnetic fields in the overall dynamics of non-rotating CCSNe. 

Our simulations further suggest that the SASI may play a role in determining the strength 
of the magnetic field in proto-neutron stars and young pulsars. We estimate that the magnetic 
energy accumulated on the PNS may account for magnetic field strengths exceeding 10^^ G. 

The evolution of magnetic fields in our simulations remains sensitive to the spatial resolution. 
The magnetic energy attained and the rate at which the magnetic energy grows increase with 
increasing grid resolution. Data extracted from our simulations suggest that the magnetic energy 
may grow exponentially on a millisecond timescale (~ Teddy) under more realistic physical 
conditions (i.e., huge magnetic Reynolds numbers), as opposed to the ~ 50-60 ms timescale 
measured directly in our runs. Model B13 attains dynamically relevant magnetic fields, which 
are also subject to significant numerical dissipation. The MHD evolution displayed by this model 
during the saturated state is probably spurious and the dynamical impact of the S-fields may 
be underestimated. We therefore caution that the uncertainties associated with the sensitivity 
to numerical resolution prevents us from completely dismissing magnetic fields as unimportant 
to the explosion dynamics of weakly rotating progenitors. 

In summary, we conclude from our simulations that magnetic fields in CCSNe may be 
amplified exponentially by turbulence driven by the spiral SASI mode. Details on the impact of 
SASI-amplified magnetic fields on explosion dynamics remain unclear, but on energetic grounds 
alone the role of these magnetic fields is likely sub-dominant. The simulations further suggest 
that small-scale PNS magnetic fields in the 10^^ — 10^^ G range may be formed, which may be 
sufficient to power some of the energetic activity that define AXPs and SGRs. 
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